Table of contents

  • Importing Libraries
  • Configuring Visualization Parameters
  • Configuring Other Notebook Parameters
  • Pre-installing Custom Functions
  • Practicing in Stages
    • DICOM in Python
      • Extraction of DICOM Files
      • Visualization of 2D DICOM Files
      • Visualization of 3D Data Stored as Multiple 2D DICOM Files
      • Extraction of DICOM Files by Auto-Detection
    • NIfTI in Python
      • Extraction of NIfTI File

Importing Libraries¶

In [1]:
# Import Pydicom before using Pydicom functions

# DICOM (Digital Imaging and Communications in Medicine) is the standard protocol for
# the management and transmission of medical images and related data, and is used by
# many healthcare facilities

# Pydicom is a pure Python package for working with DICOM files
import pydicom
In [2]:
# The `dicom2nifti` library converts DICOM files to NIfTI files and supports most anatomical CT
# and MRI data, especially MRI, which supports most 4D data such as DTI data and fMRI data
import dicom2nifti

# NiBabel is a pure Python package that supports a growing number of neuroimaging file formats,
# including the NIfTI-1 format and the NIfTI-2 format, which is an extension of the former

# NiBabel provides both format-independent access to high-level neuroimaging and format-specific APIs with
# different levels of access to all available information in a particular file format

# NiBabel's API provides full or selective access to file header information (metadata), while image data
# is provided via NumPy arrays
import nibabel as nib
In [3]:
# The `pathlib` module is similar to the `os.path` module, but `pathlib` provides a more
# advanced and convenient interface than `os.path`

# It is possible to use `pathlib` to represent file paths as specialized `Path` objects
# instead of plain strings
from pathlib import Path

# SimpleITK is a simplified programming interface to the algorithms and data structures
# of the Insight Toolkit (ITK) for segmentation, registration and advanced image analysis

# SimpleITK supports bindings for multiple programming languages including C++, Python, R,
# Java, C#, Lua, Ruby and TCL

# Combining SimpleITK’s Python bindings with the Jupyter notebook web application creates
# an environment which facilitates collaborative development of biomedical image analysis
# workflows

# In this project, SimpleITK provides the ability to automatically detect and read all DICOM
# files so that the user does not have to manage file reading or slice sorting
import SimpleITK as sitk
In [4]:
import numpy as np

import torch
from torchvision.utils import make_grid
In [5]:
import matplotlib as mpl
import matplotlib.gridspec as gridspec
import matplotlib.pyplot as plt
In [6]:
from datetime import datetime
from functools import wraps
import itertools
import math
import random
import reprlib
import sys

Configuring Visualization Parameters¶

In [7]:
%matplotlib inline
In [8]:
XINHUI = "#7a7374"
XUEBAI = "#fffef9"
YINBAI = "#f1f0ed"
YINHUI = "#918072"

figure_size = (16, 9)
In [9]:
custom_params = {
    "axes.axisbelow": True,
    "axes.edgecolor": YINBAI,
    "axes.facecolor": XUEBAI,
    "axes.grid": True,
    "axes.labelcolor": XINHUI,
    "axes.spines.right": False,
    "axes.spines.top": False,
    "axes.titlecolor": XINHUI,
    "figure.edgecolor": YINBAI,
    "figure.facecolor": XUEBAI,
    "grid.alpha": 0.8,
    "grid.color": YINBAI,
    "grid.linestyle": "--",
    "grid.linewidth": 1.2,
    "legend.edgecolor": YINHUI,
    "patch.edgecolor": XUEBAI,
    "patch.force_edgecolor": True,
    "text.color": XINHUI,
    "xtick.color": YINHUI,
    "ytick.color": YINHUI,
}

mpl.rcParams.update(custom_params)

Configuring Other Notebook Parameters¶

In [10]:
reprlib_rules = reprlib.Repr()
reprlib_rules.maxother = 250

Pre-installing Custom Functions¶

In [11]:
sys.path.append("../")
In [12]:
from Modules import *

Practicing in Stages¶

DICOM in Python¶

Extraction of DICOM Files¶

In [13]:
# The dataset used in this practice project is a very small subset of CT images extracted from
# the Cancer Imaging Archive (TCIA), which contains intermediate slices of all CT images
# in which valid age, mode, and contrast markers can be found

# The dataset is provided by the Kaggel Datasets called CT Medical Imaging
# (https://www.kaggle.com/datasets/kmader/siim-medical-images),
# license type is Attribution 3.0 Unported (CC BY 3.0)
# https://creativecommons.org/licenses/by/3.0

# The link to the TCIA archive of the full dataset is
# https://wiki.cancerimagingarchive.net/display/Public/TCGA-LUAD
dir_path = "../Datasets/Kaggle - CT Medical Images/dicom_dir/"
sample_dcm = "ID_0000_AGE_0060_CONTRAST_1_CT.dcm"
# The main function of Pydicom to read and parse DICOM files is `read_file`
dicom_file = pydicom.read_file(dir_path + sample_dcm)

tabulation = Form_Generator()
tabulation.heading_printer("Initial understanding of DICOM file")

statements = [
    """
dir_path = "../Datasets/Kaggle - CT Medical Images/dicom_dir/"
sample_dcm = "ID_0000_AGE_0060_CONTRAST_1_CT.dcm"
dicom_file = pydicom.read_file(dir_path + sample_dcm)
"""
]
tabulation.statement_generator(statements)

variables = ["dicom_file"]
values = [str(reprlib_rules.repr(dicom_file))]
tabulation.variable_generator(variables, values)

expressions = [
    "dicom_file[0x0028, 0x0010]",
    "dicom_file[0x0028, 0x0011]",
    "dicom_file[0x0018, 0x0015]",
    "dicom_file.Rows",
    "dicom_file.Columns",
    "dicom_file.BodyPartExamined",
    "dicom_file.keys()",
    "dicom_file.values()",
    "dicom_file.dir()",
    'dicom_file.dir("Image")',
]
results = [
    str(dicom_file[0x0028, 0x0010]),
    str(dicom_file[0x0028, 0x0011]),
    str(dicom_file[0x0018, 0x0015]),
    str(dicom_file.Rows),
    str(dicom_file.Columns),
    str(dicom_file.BodyPartExamined),
    str(reprlib_rules.repr(dicom_file.keys())),
    str(reprlib_rules.repr(dicom_file.values())),
    str(reprlib_rules.repr(dicom_file.dir())),
    str(reprlib_rules.repr(dicom_file.dir("Image"))),
]
tabulation.expression_generator(expressions, results, 1)
Initial understanding of DICOM file

    +-------------------------------------------------------+
    | Statement                                             |
    +-------------------------------------------------------+
    | dir_path = "../Datasets/Kaggle - CT Medical           |
    |     Images/dicom_dir/"                                |
    | sample_dcm = "ID_0000_AGE_0060_CONTRAST_1_CT.dcm"     |
    | dicom_file = pydicom.read_file(dir_path + sample_dcm) |
    +-------------------------------------------------------+
    +------------+------------------------------------------------+
    | Variable   | Value                                          |
    +------------+------------------------------------------------+
    | dicom_file | Dataset.file_meta                              |
    |            |         -------------------------------        |
    |            | (0002, 0000) File Meta Information Group       |
    |            |         Length  UL: 194                        |
    |            | (0002, 0001) Fil...Group Length                |
    |            |                 UL: 524296                     |
    |            | (7fe0, 0010) Pixel Data                        |
    |            |           OW: Array of 524288 elements         |
    +------------+------------------------------------------------+
    +-----------------------------+-------------------------------+
    | Expression                  | Result                        |
    +-----------------------------+-------------------------------+
    | dicom_file[0x0028, 0x0010]  | (0028, 0010) Rows             |
    |                             |                     US: 512   |
    | dicom_file[0x0028, 0x0011]  | (0028, 0011) Columns          |
    |                             |                     US: 512   |
    | dicom_file[0x0018, 0x0015]  | (0018, 0015) Body Part        |
    |                             |  Examined                     |
    |                             |  CS: 'CHEST'                  |
    | dicom_file.Rows             | 512                           |
    | dicom_file.Columns          | 512                           |
    | dicom_file.BodyPartExamined | CHEST                         |
    | dicom_file.keys()           | dict_keys([(0008, 0000),      |
    |                             |  (0008, 0005), (0008, 0008),  |
    |                             |  (0008, 0016), (0008, 0018),  |
    |                             |  (0008, 0020), (0008, 0021),  |
    |                             |  (0008, 0022), ...095, 0010), |
    |                             |  (0097, 0000), (0097, 0010),  |
    |                             |  (0099, 0000), (0099, 0010),  |
    |                             |  (7003, 0000), (7003, 0010),  |
    |                             |  (7fe0, 0000), (7fe0, 0010)]) |
    | dicom_file.values()         | dict_values([(0008, 0000)     |
    |                             |  Group Length                 |
    |                             |         UL: 430, (0008, 0005) |
    |                             |  Specific Character Set       |
    |                             |         CS:...up Length       |
    |                             |                   UL: 524296, |
    |                             |  (7fe0, 0010) Pixel Data      |
    |                             |                      OW:      |
    |                             |  Array of 524288 elements])   |
    | dicom_file.dir()            | ['AccessionNumber',           |
    |                             |  'AcquisitionDate',           |
    |                             |  'AcquisitionNumber',         |
    |                             |  'AcquisitionTime',           |
    |                             |  'BitsAllocated',             |
    |                             |  'BitsStored', ...]           |
    | dicom_file.dir("Image")     | ['ImageDimensions',           |
    |                             |  'ImageFormat',               |
    |                             |  'ImageGeometryType',         |
    |                             |  'ImageLocation',             |
    |                             |  'ImageOrientation',          |
    |                             |  'ImageOrientationPatient',   |
    |                             |  ...]                         |
    +-----------------------------+-------------------------------+
In [14]:
# DICOM data has an attribute called `pixel_array` that provides more useful pixel data
# for uncompressed images that can be passed to the graphics library for viewing

# To use this attribute, the system must have the NumPy numeric package installed,
# since `pixel_array` returns a NumPy array
ct = dicom_file.pixel_array

tabulation = Form_Generator()
tabulation.heading_printer("Getting pixel data from DICOM file")

statements = ["ct = dicom_file.pixel_array"]
tabulation.statement_generator(statements)

variables = ["ct"]
values = [str(reprlib_rules.repr(ct))]
tabulation.variable_generator(variables, values)

expressions = ["ct.shape"]
results = [str(ct.shape)]
tabulation.expression_generator(expressions, results)
Getting pixel data from DICOM file

    +-----------------------------+
    | Statement                   |
    +-----------------------------+
    | ct = dicom_file.pixel_array |
    +-----------------------------+
    +----------+------------------------------------------------+
    | Variable | Value                                          |
    +----------+------------------------------------------------+
    | ct       | array([[0, 0, 0, ..., 0, 0, 0],                |
    |          |        [0, 0, 0, ..., 0, 0, 0],                |
    |          |        [0, 0, 0, ..., 0, 0, 0],                |
    |          |        ...,                                    |
    |          |        [0, 0, 0, ..., 0, 0, 0],                |
    |          |        [0, 0, 0, ..., 0, 0, 0],                |
    |          |        [0, 0, 0, ..., 0, 0, 0]], dtype=uint16) |
    +----------+------------------------------------------------+
    +------------+------------+
    | Expression | Result     |
    +------------+------------+
    | ct.shape   | (512, 512) |
    +------------+------------+

Visualization of 2D DICOM Files¶

In [15]:
def image_display(image, ax, title, cmap):
    ax.imshow(image, cmap)
    ax.grid(False)
    ax.set_title(title, loc="center", pad=10)
    x_ticks = list(range(0, image.shape[1], 50))
    y_ticks = list(range(0, image.shape[0], 50))
    ax.set(xticks=x_ticks, xticklabels=x_ticks,
           yticks=y_ticks, yticklabels=y_ticks)
    ax.set_xlim(left=0)
    ax.set_ylim(top=0)
    ax.minorticks_on()
    return ax


# Path classes are divided between pure paths, which provide purely computational
# operations without I/O, and concrete paths, which inherit from pure paths but also
# provide I/O operations
path_object = Path(dir_path)
# `PurePath.name` returns a string representing the final path component, excluding
# the drive and root directory (if any)

# When a path points to a directory, `Path.iterdir` generates a path object of the directory's
# contents

# `Path.is_file` returns True if the path points to a normal file or to a symbolic link
# to a normal file, False if it points to another type of file
random_dicom_path = random.choice(
    [
        file.name
        for file in path_object.iterdir()
        if file.is_file() & (pydicom.read_file(file).BodyPartExamined != "CHEST")
    ]
)

random_dicom_file = pydicom.read_file(dir_path + random_dicom_path)

plt.rcParams["figure.figsize"] = (
    figure_size[0] / 3 * 2, figure_size[1] / 4 * 5)

fig, axs = plt.subplots(nrows=2, ncols=2)

image_display(
    ct,
    axs[0, 0],
    f"{dicom_file.BodyPartExamined.capitalize()} CT medical image displayed using "
    "colormap 'gray'",
    cmap="gray",
)

image_display(
    ct,
    axs[0, 1],
    f"{dicom_file.BodyPartExamined.capitalize()} CT medical image displayed using "
    "colormap 'bone'",
    cmap="bone",
)

image_display(
    random_dicom_file.pixel_array,
    axs[1, 0],
    f"{random_dicom_file.BodyPartExamined.capitalize()} CT medical image displayed using "
    "colormap 'gray'",
    cmap="gray",
)

image_display(
    random_dicom_file.pixel_array,
    axs[1, 1],
    f"{random_dicom_file.BodyPartExamined.capitalize()} CT medical image displayed using "
    "colormap 'bone'",
    cmap="bone",
)


fig.suptitle(
    "Visual Comparison of CT Medical Image Display Using Different Colormaps",
    fontsize="x-large",
    x=0.5,
    y=0,
)

plt.tight_layout()
plt.show()
In [16]:
def annotation_color(cmap):
    cmap_color = mpl.colormaps[cmap]
    return cmap_color(0.85)


def first_max(list):
    max_element = list[0]
    for element in list[1:]:
        if len(element) > len(max_element):
            max_element = element
    return max_element


def text_aligner(text):
    line_list = text.split("\n")
    max_length = len(first_max(line_list))
    line_list = [
        (":" + " " * (max_length - len(line))).join(line.split(":"))
        if len(line) < max_length
        else line
        for line in line_list
    ]
    text = "\n".join(line_list)
    return text


def plane_information_annotator(func):
    @wraps(func)
    def wrapper(ax, cmap, dicom_file, fontsize_adjustment):
        plane = func(ax, cmap, dicom_file, fontsize_adjustment)
        text = f"Image: {dicom_file.Columns} x {dicom_file.Rows}"
        text += f"\nImage Plane: {plane}"
        text += f"\nPatient Position: {dicom_file.PatientPosition}"
        text += f"\nSlice Thickness: {dicom_file.SliceThickness:.1f} mm"
        text += f"\nSlice Location: {dicom_file.SliceLocation:.1f} mm"
        ax.text(
            0.575,
            0.825,
            text_aligner(text),
            transform=ax.transAxes,
            fontsize=9 + fontsize_adjustment,
            fontfamily="monospace",
            c=annotation_color(cmap),
        )

    return wrapper


def orientation_annotator(func):
    @wraps(func)
    def wrapper(ax, cmap, dicom_file, fontsize_adjustment):
        plane, top, right, bottom, left = func(dicom_file)
        ax.text(
            0.4875,
            0.95,
            top,
            transform=ax.transAxes,
            fontweight="heavy",
            fontsize=11 + fontsize_adjustment,
            fontfamily="monospace",
            c=annotation_color(cmap),
        )
        ax.text(
            0.4875,
            0.025,
            bottom,
            transform=ax.transAxes,
            fontweight="heavy",
            fontsize=11 + fontsize_adjustment,
            fontfamily="monospace",
            c=annotation_color(cmap),
        )
        ax.text(
            0.95,
            0.4875,
            right,
            transform=ax.transAxes,
            fontweight="heavy",
            fontsize=11 + fontsize_adjustment,
            fontfamily="monospace",
            c=annotation_color(cmap),
        )
        ax.text(
            0.025,
            0.4875,
            left,
            transform=ax.transAxes,
            fontweight="heavy",
            fontsize=11 + fontsize_adjustment,
            fontfamily="monospace",
            c=annotation_color(cmap),
        )
        return plane

    return wrapper


def image_orientation(func):
    @wraps(func)
    def wrapper(*args, **kwargs):
        Left, Right, Anterior, Posterior, Head, Feet = "L", "R", "A", "P", "H", "F"
        axial_plane = [Anterior, Right, Posterior, Left]
        coronal_plane = [Head, Left, Feet, Right]
        sagittal_plane = [Head, Anterior, Feet, Posterior]
        plane = func(*args, **kwargs)
        # Patient Position (0018,5100) specifies the position of the patient relative to
        # the space of the imaging device and is used for annotation purposes only,
        # but does not provide a precise mathematical relationship between the patient and
        # the imaging device
        patient_position = dicom_file.PatientPosition
        if plane == "Axial":
            term_list = ["S", "DR", "P", "DL"]
            top = 4 - term_list.index(patient_position[2:])
            if patient_position.startswith("HF"):
                right = top + 1
                left = right + 2
            else:
                left = top + 1
                right = left + 2
            bottom = top + 2
            top, right, bottom, left = [
                axial_plane[index] if index < 4 else axial_plane[index - 4]
                for index in [top, right, bottom, left]
            ]
        elif plane == "Coronal":
            term_list = ["HF", "RF", "FF", "LF"]
            top = 4 - term_list.index(patient_position[:2])
            if patient_position.endswith("S"):
                right = top + 1
                left = right + 2
            else:
                left = top + 1
                right = left + 2
            bottom = top + 2
            top, right, bottom, left = [
                coronal_plane[index] if index < 4 else coronal_plane[index - 4]
                for index in [top, right, bottom, left]
            ]
        elif plane == "Sagittal":
            term_list = ["HF", "PF", "FF", "AF"]
            top = 4 - term_list.index(patient_position[:2])
            if patient_position.endswith("DL"):
                right = top + 1
                left = right + 2
            else:
                left = top + 1
                right = left + 2
            bottom = top + 2
            top, right, bottom, left = [
                sagittal_plane[index] if index < 4 else sagittal_plane[index - 4]
                for index in [top, right, bottom, left]
            ]
        else:
            top, right, bottom, left = np.repeat("", 4)
        return plane, top, right, bottom, left

    return wrapper


@plane_information_annotator
@orientation_annotator
@image_orientation
def get_plane_information(dicom_file):
    # The DICOM standard defines a patient-oriented reference coordinate system (RCS) that
    # enables the user to measure the position and orientation of the image relative to
    # the patient

    # Image Orientation (Patient) (0020,0037) specifies the cosine of the orientation of
    # the first row and column relative to the patient, which should be provided as a pair,
    # where the row values for the x, y, and z axes are followed by the column values for
    # the x, y, and z axes

    # The orientation of the axes is determined solely by the orientation of the patient, i.e.,
    # the RCS and Image Orientation (Patient) (0020,0037) values specify the orientation of
    # the image frame rows and columns
    IOP = dicom_file.ImageOrientationPatient
    row_xyz = IOP[:3]
    column_xyz = IOP[3:]
    # `numpy.cross` returns the cross product of two vector arrays

    # The geometric interpretation of the cross product is that two vectors determine a plane,
    # and the cross product points in a direction different from both vectors

    # The Patient-Based Coordinate System is a right-handed system, i.e., the vector
    # cross product of a unit vector along the positive x-axis and a unit vector along
    # the positive y-axis is equal to a unit vector along the positive z-axis
    plane = [abs(round(x)) for x in np.cross(row_xyz, column_xyz)]
    # The Image Type (0008,0008) identifies important image identification features
    # and it returns the Pixel Data Characteristics, Patient Examination Characteristics,
    # Modality Specific Characteristics, and Implementation Specific Identifiers

    # While the third value (Modality Specific Characteristics) should identify any
    # specialization specific to the Image Information Object Definition (IOD), it
    # and the values that follow it are not mandatory, unlike the first two values,
    # and therefore some DICOM data does not have this value
    if plane[0] == 1 and plane[1] == 0 and plane[2] == 0:
        return "Sagittal"
    elif plane[0] == 0 and plane[1] == 1 and plane[2] == 0:
        return "Coronal"
    elif plane[0] == 0 and plane[1] == 0 and plane[2] == 1:
        return "Axial"
    else:
        "Unknown"


def get_patient_information(ax, cmap, dicom_file, fontsize_adjustment, patien_age):
    text = f"Patient ID: {dicom_file.PatientID}"
    text += f"\nPatient's Sex: {dicom_file.PatientSex}"
    try:
        patien_age = dicom_file.PatientAge.strip(str(0))[:-1]
    except AttributeError:
        patien_age = patien_age
    text += f"\nPatient's Age: {patien_age} y"
    text += f"\nModality: {dicom_file.Modality}"
    text += f"\nBody Part Examined: {dicom_file.BodyPartExamined}"
    ax.text(
        0.04,
        0.825,
        text_aligner(text),
        transform=ax.transAxes,
        fontsize=9 + fontsize_adjustment,
        fontfamily="monospace",
        c=annotation_color(cmap),
    )


def get_date_information(ax, cmap, dicom_file, fontsize_adjustment):
    date = datetime.strptime(dicom_file.StudyDate, "%Y%m%d").strftime("%Y.%m.%d")
    text = f"Study Date: {date}"
    ax.text(
        0.04,
        0.05,
        text,
        transform=ax.transAxes,
        fontsize=7 + fontsize_adjustment,
        fontfamily="monospace",
        c=annotation_color(cmap),
    )


def get_manufacturer_information(ax, cmap, dicom_file, fontsize_adjustment):
    text = f"Manufacturer: {dicom_file.Manufacturer}"
    ax.text(
        0.525,
        0.05,
        text.rjust(36),
        transform=ax.transAxes,
        fontsize=7 + fontsize_adjustment,
        fontfamily="monospace",
        c=annotation_color(cmap),
    )


def DICOM_2D_image(ax, cmap, dicom_file, title, fontsize_adjustment=0, patien_age="__"):
    ax.imshow(dicom_file.pixel_array, cmap)
    ax.grid(False)
    ax.set_title(title, loc="center", pad=10, fontsize=12 + fontsize_adjustment)
    ax.set(xticks=[], yticks=[], frame_on=False)
    get_patient_information(ax, cmap, dicom_file, fontsize_adjustment, patien_age)
    get_plane_information(ax, cmap, dicom_file, fontsize_adjustment)
    get_date_information(ax, cmap, dicom_file, fontsize_adjustment)
    get_manufacturer_information(ax, cmap, dicom_file, fontsize_adjustment)


fig, axs = plt.subplots(2, 2, figsize=(figure_size[0] / 3 * 2, figure_size[1] / 4 * 5))

DICOM_2D_image(
    axs[0, 0],
    "gray",
    dicom_file,
    f"{dicom_file.BodyPartExamined.capitalize()} CT medical image displayed using "
    "colormap 'gray'",
)

DICOM_2D_image(
    axs[0, 1],
    "bone",
    dicom_file,
    f"{dicom_file.BodyPartExamined.capitalize()} CT medical image displayed using "
    "colormap 'bone'",
)

DICOM_2D_image(
    axs[1, 0],
    "gray",
    random_dicom_file,
    f"{random_dicom_file.BodyPartExamined.capitalize()} CT medical image displayed using "
    "colormap 'gray'",
)

DICOM_2D_image(
    axs[1, 1],
    "bone",
    random_dicom_file,
    f"{random_dicom_file.BodyPartExamined.capitalize()} CT medical image displayed using "
    "colormap 'bone'",
)


fig.suptitle(
    "Visualization of CT Medical Image Display with Supplementary Information",
    fontsize="x-large",
    x=0.5,
    y=0,
)

plt.tight_layout()
plt.show()

Visualization of 3D Data Stored as Multiple 2D DICOM Files¶

In [17]:
# The dataset used in this practical project is the MRI DICOM dataset of the head of
# a 52-year-old normal male mathematics professor

# The subject of the experiment is the author, who suffers from small vertical strabismus
# (hypertropia), which is a misalignment of the eyes, which can be seen in this dataset

# This dataset was provided by Lionheart, William R.B. in 2015, available online
# (https://zenodo.org/record/16956 or http://doi.org/10.5281/zenodo.16956),
# license type is Attribution-ShareAlike 4.0 International (CC BY-SA 4.0)
# https://creativecommons.org/licenses/by-sa/4.0
path_to_head_mri = Path(
    "../Datasets/An MRI Dicom Data Set of the Head of a Normal Male Human Aged 52/"
    "ST000000/SE000001/"
)
# `Path.glob` selects the given relative pattern in the directory represented by the given path
# and yields all matching files (of any type)

# Since `Path.glob` returns a generator, it is converted to a list here for easy display
all_files = list(path_to_head_mri.glob("*"))

mri_data = []

for path in all_files:
    # Read the individual DICOM files one by one in the list returned by the path generator
    data = pydicom.read_file(path)
    mri_data.append(data)

tabulation = Form_Generator()
tabulation.heading_printer(
    "Reading all DICOM files that match the specified pattern")

statements = [
    """
path_to_head_mri = Path(
    "../Datasets/An MRI Dicom Data Set of the Head of a Normal Male Human Aged 52/"
    "ST000000/SE000001/"
)
all_files = list(path_to_head_mri.glob("*"))

mri_data = []

for path in all_files:
    data = pydicom.read_file(path)
    mri_data.append(data)
"""
]
tabulation.statement_generator(statements)

variables = ["str(path_to_head_mri)", "all_files"]
values = [str(path_to_head_mri), str(reprlib_rules.repr(all_files))]
tabulation.variable_generator(variables, values, 1)

expressions = [
    "all_files[0].parts",
    "all_files[0].parent",
    "all_files[0].name",
    "all_files[0].suffixes",
    "all_files[0].stem",
    "len(all_files)",
    "mri_data[0]",
    "len(mri_data)",
]
results = [
    str(all_files[0].parts),
    str(all_files[0].parent),
    str(all_files[0].name),
    str(all_files[0].suffixes),
    str(all_files[0].stem),
    str(len(all_files)),
    str(reprlib_rules.repr(mri_data[0])),
    str(len(mri_data)),
]
tabulation.expression_generator(expressions, results, 1)
Reading all DICOM files that match the specified pattern

    +---------------------------------------------------------+
    | Statement                                               |
    +---------------------------------------------------------+
    | path_to_head_mri = Path(                                |
    |     "../Datasets/An MRI Dicom Data Set of the Head of a |
    |     Normal Male Human Aged 52/"                         |
    |     "ST000000/SE000001/"                                |
    | )                                                       |
    | all_files = list(path_to_head_mri.glob("*"))            |
    |                                                         |
    | mri_data = []                                           |
    |                                                         |
    | for path in all_files:                                  |
    |     data = pydicom.read_file(path)                      |
    |     mri_data.append(data)                               |
    +---------------------------------------------------------+
    +-----------------------+-------------------------------------+
    | Variable              | Value                               |
    +-----------------------+-------------------------------------+
    | str(path_to_head_mri) | ../Datasets/An MRI Dicom Data Set   |
    |                       |  of the Head of a Normal Male Human |
    |                       |  Aged 52/ST000000/SE000001          |
    | all_files             | [PosixPath('../Datasets/An MRI      |
    |                       |  Dicom Data Set of the Head of a    |
    |                       |  Normal Male Human Aged             |
    |                       |  52/ST000000/SE000001/MR000000'),   |
    |                       |  PosixPath('../Datasets/An MRI      |
    |                       |  Dicom Data Set of the Head of a    |
    |                       |  Normal Male Human Aged             |
    |                       |  52/ST000000/SE000001/MR000007'),   |
    |                       |  PosixPath('../Datasets/An MRI      |
    |                       |  Dicom Data Set of the Head of a    |
    |                       |  Normal Male Human Aged             |
    |                       |  52/ST000000/SE000001/MR000009'),   |
    |                       |  PosixPath('../Datasets/An MRI      |
    |                       |  Dicom Data Set of the Head of a    |
    |                       |  Normal Male Human Aged             |
    |                       |  52/ST000000/SE000001/MR000008'),   |
    |                       |  PosixPath('../Datasets/An MRI      |
    |                       |  Dicom Data Set of the Head of a    |
    |                       |  Normal Male Human Aged             |
    |                       |  52/ST000000/SE000001/MR000006'),   |
    |                       |  PosixPath('../Datasets/An MRI      |
    |                       |  Dicom Data Set of the Head of a    |
    |                       |  Normal Male Human Aged             |
    |                       |  52/ST000000/SE000001/MR000001'),   |
    |                       |  ...]                               |
    +-----------------------+-------------------------------------+
    +-----------------------+-------------------------------------+
    | Expression            | Result                              |
    +-----------------------+-------------------------------------+
    | all_files[0].parts    | ('..', 'Datasets', 'An MRI Dicom    |
    |                       |  Data Set of the Head of a Normal   |
    |                       |  Male Human Aged 52', 'ST000000',   |
    |                       |  'SE000001', 'MR000000')            |
    | all_files[0].parent   | ../Datasets/An MRI Dicom Data Set   |
    |                       |  of the Head of a Normal Male Human |
    |                       |  Aged 52/ST000000/SE000001          |
    | all_files[0].name     | MR000000                            |
    | all_files[0].suffixes | []                                  |
    | all_files[0].stem     | MR000000                            |
    | len(all_files)        | 27                                  |
    | mri_data[0]           | Dataset.file_meta                   |
    |                       |  -------------------------------    |
    |                       | (0002, 0000) File Meta Information  |
    |                       |  Group Length  UL: 214              |
    |                       | (0002, 0001) Fil...Group Length     |
    |                       |                     UL: 131084      |
    |                       | (7fe0, 0010) Pixel Data             |
    |                       |               OW: Array of 131072   |
    |                       |  elements                           |
    | len(mri_data)         | 27                                  |
    +-----------------------+-------------------------------------+
In [18]:
# DICOM files may not be sorted according to their actual image positions, this can
# be verified by checking Slice Location (0020,1041)

# Slice Location (0020,1041) identifies the relative position of the image plane in millimeters
unordered_slices = {
    index: float(slice.SliceLocation) for index, slice in enumerate(mri_data)
}

# For better viewing and processing of 3D data stored with multiple 2D DICOM files,
# these files must be sorted; otherwise the complete scanned image will become cluttered
# and useless
mri_data_ordered = sorted(mri_data, key=lambda slice: slice.SliceLocation)

ordered_slices = {
    index: float(slice.SliceLocation) for index, slice in enumerate(mri_data_ordered)
}


tabulation = Form_Generator()
tabulation.heading_printer("Sorting of read DICOM files")

statements = [
    """
unordered_slices = {
    index: float(slice.SliceLocation) for index, slice in enumerate(mri_data)
}

mri_data_ordered = sorted(mri_data, key=lambda slice: slice.SliceLocation)

ordered_slices = {
    index: float(slice.SliceLocation) for index, slice in enumerate(mri_data_ordered)
}
"""
]
tabulation.statement_generator(statements)

variables = ["unordered_slices", "ordered_slices"]
values = [
    str(unordered_slices),
    str(ordered_slices),
]
tabulation.variable_generator(variables, values, 1)

expressions = [
    "mri_data_ordered[0]",
    "len(mri_data_ordered)",
    "len(unordered_slices)",
    "len(ordered_slices)",
]
results = [
    str(reprlib_rules.repr(mri_data_ordered[0])),
    str(len(mri_data_ordered)),
    str(len(unordered_slices)),
    str(len(ordered_slices)),
]
tabulation.expression_generator(expressions, results, 1)
Sorting of read DICOM files

    +-----------------------------------------------------------+
    | Statement                                                 |
    +-----------------------------------------------------------+
    | unordered_slices = {                                      |
    |     index: float(slice.SliceLocation) for index, slice in |
    |     enumerate(mri_data)                                   |
    | }                                                         |
    |                                                           |
    | mri_data_ordered = sorted(mri_data, key=lambda slice:     |
    |     slice.SliceLocation)                                  |
    |                                                           |
    | ordered_slices = {                                        |
    |     index: float(slice.SliceLocation) for index, slice in |
    |     enumerate(mri_data_ordered)                           |
    | }                                                         |
    +-----------------------------------------------------------+
    +------------------+------------------------------------------+
    | Variable         | Value                                    |
    +------------------+------------------------------------------+
    | unordered_slices | {0: 0.0, 1: 41.9999963629367, 2:         |
    |                  |  53.9999958207213, 3: 47.9999970362677,  |
    |                  |  4: 35.9999959546749, 5:                 |
    |                  |  5.99999663091323, 6: 137.999998321624,  |
    |                  |  7: 143.999998928727, 8:                 |
    |                  |  71.9999961590453, 9: 89.9999955528687,  |
    |                  |  10: 83.9999967682912, 11:               |
    |                  |  77.999996227574, 12: 149.999999502083,  |
    |                  |  13: 131.999997780749, 14:               |
    |                  |  23.9999946081714, 15: 17.9999979772582, |
    |                  |  16: 11.9999973042441, 17:               |
    |                  |  29.9999952815023, 18: 107.999995419197, |
    |                  |  19: 119.999996566542, 20:               |
    |                  |  95.9999960937442, 21: 65.9999961939969, |
    |                  |  22: 59.9999962290673, 23:               |
    |                  |  101.999994745866, 24: 125.999997173645, |
    |                  |  25: 155.999992554172, 26:               |
    |                  |  113.999995959439}                       |
    | ordered_slices   | {0: 0.0, 1: 5.99999663091323, 2:         |
    |                  |  11.9999973042441, 3: 17.9999979772582,  |
    |                  |  4: 23.9999946081714, 5:                 |
    |                  |  29.9999952815023, 6: 35.9999959546749,  |
    |                  |  7: 41.9999963629367, 8:                 |
    |                  |  47.9999970362677, 9: 53.9999958207213,  |
    |                  |  10: 59.9999962290673, 11:               |
    |                  |  65.9999961939969, 12: 71.9999961590453, |
    |                  |  13: 77.999996227574, 14:                |
    |                  |  83.9999967682912, 15: 89.9999955528687, |
    |                  |  16: 95.9999960937442, 17:               |
    |                  |  101.999994745866, 18: 107.999995419197, |
    |                  |  19: 113.999995959439, 20:               |
    |                  |  119.999996566542, 21: 125.999997173645, |
    |                  |  22: 131.999997780749, 23:               |
    |                  |  137.999998321624, 24: 143.999998928727, |
    |                  |  25: 149.999999502083, 26:               |
    |                  |  155.999992554172}                       |
    +------------------+------------------------------------------+
    +-----------------------+-------------------------------------+
    | Expression            | Result                              |
    +-----------------------+-------------------------------------+
    | mri_data_ordered[0]   | Dataset.file_meta                   |
    |                       |  -------------------------------    |
    |                       | (0002, 0000) File Meta Information  |
    |                       |  Group Length  UL: 214              |
    |                       | (0002, 0001) Fil...Group Length     |
    |                       |                     UL: 131084      |
    |                       | (7fe0, 0010) Pixel Data             |
    |                       |               OW: Array of 131072   |
    |                       |  elements                           |
    | len(mri_data_ordered) | 27                                  |
    | len(unordered_slices) | 27                                  |
    | len(ordered_slices)   | 27                                  |
    +-----------------------+-------------------------------------+
In [19]:
# Fill the 3D array in a slice-per-slice manner
full_volume = [slice.pixel_array for slice in mri_data_ordered]
full_volume = np.array(full_volume)

tabulation = Form_Generator()
tabulation.heading_printer(
    "One-time extraction of the overall 3D pixel array for all slices"
)

statements = [
    """
full_volume = [slice.pixel_array for slice in mri_data_ordered]
full_volume = np.array(full_volume)
"""
]
tabulation.statement_generator(statements)

variables = ["full_volume"]
values = [str(reprlib_rules.repr(full_volume))]
tabulation.variable_generator(variables, values)

expressions = ["full_volume.shape", "full_volume.dtype"]
results = [str(full_volume.shape), str(full_volume.dtype)]
tabulation.expression_generator(expressions, results)
One-time extraction of the overall 3D pixel array for all slices

    +-----------------------------------------------+
    | Statement                                     |
    +-----------------------------------------------+
    | full_volume = [slice.pixel_array for slice in |
    |     mri_data_ordered]                         |
    | full_volume = np.array(full_volume)           |
    +-----------------------------------------------+
    +-------------+------------------------------------+
    | Variable    | Value                              |
    +-------------+------------------------------------+
    | full_volume | array([[[0, 0, 0, ..., 0, 0, 0],   |
    |             |         [0, 0, 0, ..., 0, 0, 0],   |
    |             |         [0, 0, 0, ..., 0, 0, 0],   |
    |             |         ...,                       |
    |             |         [0,...     ...,            |
    |             |         [0, 0, 0, ..., 0, 0, 0],   |
    |             |         [0, 0, 0, ..., 0, 0, 0],   |
    |             |         [0, 0, 0, ..., 0, 0, 0]]], |
    |             |         dtype=uint16)              |
    +-------------+------------------------------------+
    +-------------------+----------------+
    | Expression        | Result         |
    +-------------------+----------------+
    | full_volume.shape | (27, 256, 256) |
    | full_volume.dtype | uint16         |
    +-------------------+----------------+
In [20]:
def sort_checker(list_1):
    list_2 = list_1[1:]
    if all(a <= b for a, b in zip(list_1, list_2)):
        return "ascending"
    elif all(a >= b for a, b in zip(list_1, list_2)):
        return "descending"
    elif all(a == b for a, b in zip(list_1, list_2)):
        return "identical"
    else:
        return "unordered"


def title_adder(func):
    @wraps(func)
    def wrapper(inputs, ax, dict_slices, title=None, **kwargs):
        is_sorted = sort_checker(list(dict_slices.values()))
        if title is None:
            title = f"DICOM image grid with {is_sorted} slice locations"
        func(inputs, ax, dict_slices, **kwargs)
        ax.set_title(
            title,
            loc="center",
            pad=10,
        )

    return wrapper


def grid_annotator(func):
    @wraps(func)
    def wrapper(inputs, ax, dict_slices, fontsize_adjustment=0, **kwargs):
        x_positions_1, x_positions_2, y_positions, cmap = func(inputs, ax, **kwargs)
        for (y, x_1), key in zip(
            itertools.product(y_positions, x_positions_1), dict_slices.keys()
        ):
            text = f"No.{(key + 1):>03}"
            ax.text(
                x_1,
                y,
                text,
                transform=ax.transData,
                fontsize=7 + fontsize_adjustment,
                fontfamily="monospace",
                c=annotation_color(cmap),
            )
        for (y, x_2), value in zip(
            itertools.product(y_positions, x_positions_2), dict_slices.values()
        ):
            text = f"{value:>5.1f} mm"
            ax.text(
                x_2,
                y,
                text,
                transform=ax.transData,
                fontsize=7 + fontsize_adjustment,
                fontfamily="monospace",
                c=annotation_color(cmap),
            )

    return wrapper


@title_adder
@grid_annotator
def grid_DICOM_2D_image(inputs, ax, n=None, row_size=9, pad_value=1):
    if n is None:
        n = inputs.shape[0]
    nrows = math.ceil(n / row_size)
    cmap = "bone"

    # The only types supported by `torch.Tensor` are float64, float32, float16, complex64,
    # complex128, int64, int32, int16, int8, uint8 and bool, not uint16
    inputs = torch.Tensor(inputs.astype(int)).unsqueeze(1)

    x_positions_1 = [
        math.ceil(
            x * inputs[:row_size].shape[3]
            + 0.075 * inputs[:row_size].shape[3]
            + pad_value * x
        )
        for x in range(inputs[:row_size].shape[0])
    ]
    x_positions_2 = [
        math.ceil(
            x * inputs[:row_size].shape[3]
            + 0.55 * inputs[:row_size].shape[3]
            + pad_value * x
        )
        for x in range(inputs[:row_size].shape[0])
    ]
    y_positions = [
        math.ceil(
            y * inputs[:row_size].shape[2]
            + 0.125 * inputs[:row_size].shape[2]
            + pad_value * y
        )
        for y in range(nrows)
    ]

    for i in range(nrows):
        if len(inputs) > row_size:
            row_images = inputs[:row_size]
            inputs = inputs[row_size:]
        else:
            row_images = inputs[:]
        # Unlike the images in the MNIST dataset, the DICOM file images do not take values
        # in the range (0, 1), so the parameter `normalize` needs to be set to True in order to
        # move their values proportionally to the range (0, 1)

        # By querying the source code of `vision/torchvision/utils.py`, it is clear that any
        # single-channel image is converted to a three-channel image, and all are represented
        # as tensors, as shown below:
        # > if tensor.dim() == 2:  # single image H x W
        # >     tensor = tensor.unsqueeze(0)
        # > if tensor.dim() == 3:  # single image
        # >     if tensor.size(0) == 1:  # if single-channel, convert to 3-channel
        # >         tensor = torch.cat((tensor, tensor, tensor), 0)
        # >     tensor = tensor.unsqueeze(0)
        # >
        # > if tensor.dim() == 4 and tensor.size(1) == 1:  # single-channel images
        # >     tensor = torch.cat((tensor, tensor, tensor), 1)
        grid_row = make_grid(
            row_images, nrow=row_size, normalize=True, pad_value=pad_value
        )
        scalar_row = grid_row.numpy().transpose(1, 2, 0)[:, :, 0]
        if i == 0:
            if nrows != 1:
                scalar_image = scalar_row
            else:
                scalar_image = np.full(
                    (
                        scalar_row.shape[0],
                        math.ceil(scalar_row.shape[1] / inputs.shape[0] * row_size),
                    ),
                    np.nan,
                )
                scalar_image[:, : scalar_row.shape[1]] = scalar_row
        elif scalar_image.shape[1] == scalar_row.shape[1]:
            scalar_image = np.concatenate((scalar_image, scalar_row), axis=0)
        else:
            image_to_concat = np.full(
                (scalar_image.shape[0] + scalar_row.shape[0], scalar_image.shape[1]),
                np.nan,
            )
            image_to_concat[: scalar_image.shape[0], :] = scalar_image
            image_to_concat[scalar_image.shape[0] :, : scalar_row.shape[1]] = scalar_row
            scalar_image = image_to_concat
    ax.imshow(scalar_image, cmap=cmap)
    ax.set(xticks=[], yticks=[], frame_on=False)
    ax.grid(False)
    return x_positions_1, x_positions_2, y_positions, cmap


unordered_full_volume = [slice.pixel_array for slice in mri_data]
unordered_full_volume = np.array(unordered_full_volume)
descending_slices = {
    index: float(slice.SliceLocation)
    for index, slice in enumerate(mri_data_ordered[::-1])
}

fig, axs = plt.subplots(3, 1, figsize=(figure_size[0], figure_size[1] / 7 * 9))
grid_DICOM_2D_image(unordered_full_volume, axs[0], unordered_slices)
grid_DICOM_2D_image(full_volume, axs[1], ordered_slices)
grid_DICOM_2D_image(full_volume[::-1], axs[2], descending_slices)

fig.suptitle(
    "Visual Comparison of DICOM Image Grids with Unordered and Ordered Slice Locations",
    fontsize="x-large",
    x=0.5,
    y=0,
)

plt.tight_layout()
plt.show()
In [21]:
fig, axs = plt.subplots(9, 3, figsize=(
    figure_size[0] / 3 * 2, figure_size[1] / 2 * 7))

for slice_counter, (i, j) in enumerate(itertools.product(range(9), range(3))):
    DICOM_2D_image(
        axs[i][j],
        "bone",
        mri_data_ordered[slice_counter],
        f"DICOM image with {mri_data_ordered[slice_counter].SliceLocation:^.1f} mm "
        "slice location",
        fontsize_adjustment=-3,
        patien_age="52",
    )

fig.suptitle(
    "Visualization of MRI Medical Images with Information Displayed in Order of Slice Location",
    fontsize="x-large",
    x=0.5,
    y=0,
)

plt.tight_layout()
plt.show()

Extraction of DICOM Files by Auto-Detection¶

In [22]:
# The path object must be converted to a string in order for the SimpleITK package to process
# the path information
path_string = str(path_to_head_mri)
# For some image formats such as DICOM, the images also contain associated metadata,
# which is not loaded by the reader by default for time-saving reasons

# `ImageSeriesReader` class provides the ability to read a series of image files into a
# SimpleITK image, and once a series of images has been read, the metadata can be accessed
# directly from the reader

# `GetGDCMSeriesIDs` provides the ability to get all seriesIDs from a DICOM dataset
series_ids = sitk.ImageSeriesReader.GetGDCMSeriesIDs(path_string)
# `GetGDCMSeriesFileNames` generates a series of file names from the directory containing
# the DICOM dataset and series ID

# Note that the resulting series of file names (which are actually file paths) are
# automatically sorted
series_file_names = sitk.ImageSeriesReader.GetGDCMSeriesFileNames(
    path_string, series_ids[0]
)

series_reader = sitk.ImageSeriesReader()
series_reader.SetFileNames(series_file_names)

# Finally, the reader is executed by calling `Execute` to get the required data
image_data = series_reader.Execute()

tabulation = Form_Generator()
tabulation.heading_printer(
    "Automatic recognition and recall of series file names for DICOM images"
)

statements = [
    """
path_string = str(path_to_head_mri)
series_ids = sitk.ImageSeriesReader.GetGDCMSeriesIDs(path_string)
series_file_names = sitk.ImageSeriesReader.GetGDCMSeriesFileNames(
    path_string, series_ids[0]
)

series_reader = sitk.ImageSeriesReader()
series_reader.SetFileNames(series_file_names)

image_data = series_reader.Execute()
"""
]
tabulation.statement_generator(statements)

variables = ["series_ids", "series_file_names", "series_reader", "image_data"]
values = [
    str(reprlib_rules.repr(series_ids)),
    str(reprlib_rules.repr(series_file_names)),
    str(reprlib_rules.repr(series_reader)),
    str(reprlib_rules.repr(image_data)),
]
tabulation.variable_generator(variables, values, 1)

expressions = [
    "type(series_ids)",
    "len(series_ids)",
    "type(series_file_names)",
    "len(series_file_names)",
    "len(image_data)",
    "image_data.GetSize()",
]
results = [
    str(type(series_ids)),
    str(len(series_ids)),
    str(type(series_file_names)),
    str(len(series_file_names)),
    str(len(image_data)),
    str(image_data.GetSize()),
]
tabulation.expression_generator(expressions, results)
Automatic recognition and recall of series file names for DICOM images

    +----------------------------------------------------------+
    | Statement                                                |
    +----------------------------------------------------------+
    | path_string = str(path_to_head_mri)                      |
    | series_ids =                                             |
    |     sitk.ImageSeriesReader.GetGDCMSeriesIDs(path_string) |
    | series_file_names =                                      |
    |     sitk.ImageSeriesReader.GetGDCMSeriesFileNames(       |
    |     path_string, series_ids[0]                           |
    | )                                                        |
    |                                                          |
    | series_reader = sitk.ImageSeriesReader()                 |
    | series_reader.SetFileNames(series_file_names)            |
    |                                                          |
    | image_data = series_reader.Execute()                     |
    +----------------------------------------------------------+
    +-------------------+-----------------------------------------+
    | Variable          | Value                                   |
    +-------------------+-----------------------------------------+
    | series_ids        | ('1.3.46.67058...1413262801702',)       |
    | series_file_names | ('../Datasets/...0001/MR000000',        |
    |                   |  '../Datasets/...0001/MR000001',        |
    |                   |  '../Datasets/...0001/MR000002',        |
    |                   |  '../Datasets/...0001/MR000003',        |
    |                   |  '../Datasets/...0001/MR000004',        |
    |                   |  '../Datasets/...0001/MR000005', ...)   |
    | series_reader     | <SimpleITK.SimpleITK.ImageSeriesReader; |
    |                   |  proxy of <Swig Object of type          |
    |                   |  'itk::simple::ImageSeriesReader *' at  |
    |                   |  0x2855a3db0> >                         |
    | image_data        | <SimpleITK.SimpleITK.Image; proxy of    |
    |                   |  <Swig Object of type 'std::vector<     |
    |                   |  itk::simple::Image >::value_type *' at |
    |                   |  0x294029d70> >                         |
    +-------------------+-----------------------------------------+
    +-------------------------+-----------------+
    | Expression              | Result          |
    +-------------------------+-----------------+
    | type(series_ids)        | ⟨class 'tuple'⟩ |
    | len(series_ids)         | 1               |
    | type(series_file_names) | ⟨class 'tuple'⟩ |
    | len(series_file_names)  | 27              |
    | len(image_data)         | 1769472         |
    | image_data.GetSize()    | (256, 256, 27)  |
    +-------------------------+-----------------+
In [23]:
# As shown above, the shape of the original SimpleITK image is (256, 256, 27) instead of
# the common shape (27, 256, 256), this is due to the different order of the image sizes,
# so it's necessary to perform the last step to convert the SimpleITK image to a NumPy array

# `GetArrayFromImage` returns a copy of the image data, which can be freely modified
# as it has no effect on the original SimpleITK image
head_mri = sitk.GetArrayFromImage(image_data)

tabulation = Form_Generator()
tabulation.heading_printer("Getting the NumPy array from the SimpleITK image")

statements = ["head_mri = sitk.GetArrayFromImage(image_data)"]
tabulation.statement_generator(statements)

variables = ["head_mri"]
values = [str(reprlib_rules.repr(head_mri))]
tabulation.variable_generator(variables, values)

expressions = [
    "type(head_mri)",
    "len(head_mri)",
    "head_mri.shape",
]
results = [
    str(type(head_mri)),
    str(len(head_mri)),
    str(head_mri.shape),
]
tabulation.expression_generator(expressions, results)
Getting the NumPy array from the SimpleITK image

    +-----------------------------------------------+
    | Statement                                     |
    +-----------------------------------------------+
    | head_mri = sitk.GetArrayFromImage(image_data) |
    +-----------------------------------------------+
    +----------+--------------------------------------------------+
    | Variable | Value                                            |
    +----------+--------------------------------------------------+
    | head_mri | array([[[0, 0, 0, ..., 0, 0, 0],                 |
    |          |         [0, 0, 0, ..., 0, 0, 0],                 |
    |          |         [0, 0, 0, ..., 0, 0, 0],                 |
    |          |         ...,                                     |
    |          |         [0,...     ...,                          |
    |          |         [0, 0, 0, ..., 0, 0, 0],                 |
    |          |         [0, 0, 0, ..., 0, 0, 0],                 |
    |          |         [0, 0, 0, ..., 0, 0, 0]]], dtype=uint16) |
    +----------+--------------------------------------------------+
    +----------------+-------------------------+
    | Expression     | Result                  |
    +----------------+-------------------------+
    | type(head_mri) | ⟨class 'numpy.ndarray'⟩ |
    | len(head_mri)  | 27                      |
    | head_mri.shape | (27, 256, 256)          |
    +----------------+-------------------------+
In [24]:
head_mri_directory_path = Path(
    "../Datasets/An MRI Dicom Data Set of the Head of a Normal Male Human Aged 52/"
    "ST000000/"
)
all_paths = sorted(list(head_mri_directory_path.glob("SE*")))
path = all_paths[-1]
new_series_ids = sitk.ImageSeriesReader.GetGDCMSeriesIDs(str(path))
new_series_file_names = sitk.ImageSeriesReader.GetGDCMSeriesFileNames(
    str(path), new_series_ids[0]
)
new_series_reader = sitk.ImageSeriesReader()
new_series_reader.SetFileNames(new_series_file_names)
new_image_data = new_series_reader.Execute()
new_head_mri = sitk.GetArrayFromImage(new_image_data)
# It is easy to see that although the ordered image data can be easily obtained by SimpleITK,
# if other related metadata is needed in addition to the image data, it is still necessary
# to obtain it in the conventional manner
new_all_files = list(path.glob("*"))
new_mri_data = sorted(
    [pydicom.read_file(p) for p in new_all_files],
    key=lambda slice: slice.SliceLocation,
)
new_ordered_slices = {
    index: float(slice.SliceLocation) for index, slice in enumerate(new_mri_data)
}

fig = plt.figure(figsize=(figure_size[0], figure_size[1]), constrained_layout=True)

gs = gridspec.GridSpec(
    nrows=2, ncols=1, figure=fig, wspace=0.08, hspace=0.04, height_ratios=[2, 3]
)

ax_1 = fig.add_subplot(gs[0, 0])
grid_DICOM_2D_image(
    head_mri, ax_1, ordered_slices, title="DICOM image grid for MRI axial plane"
)
ax_2 = fig.add_subplot(gs[1, 0])
grid_DICOM_2D_image(
    new_head_mri,
    ax_2,
    new_ordered_slices,
    title="DICOM image grid for MRI coronal plane",
    row_size=8,
    fontsize_adjustment=1,
)

fig.suptitle(
    "Visualization Comparison of Different Plane DICOM Image Grids",
    fontsize="x-large",
    x=0.5,
    y=-0.02,
)

plt.show()
In [25]:
# Create a new function based on the one used earlier, just to get the independent
# plane information
def get_plane(dicom_file):
    IOP = dicom_file.ImageOrientationPatient
    row_xyz = IOP[:3]
    column_xyz = IOP[3:]
    plane = [abs(round(x)) for x in np.cross(row_xyz, column_xyz)]
    if plane[0] == 1 and plane[1] == 0 and plane[2] == 0:
        return "sagittal"
    elif plane[0] == 0 and plane[1] == 1 and plane[2] == 0:
        return "coronal"
    elif plane[0] == 0 and plane[1] == 0 and plane[2] == 1:
        return "axial"
    else:
        "Unknown"


path = all_paths[0]
all_files_0 = list(path.glob("*"))
mri_data_0 = sorted(
    [pydicom.read_file(p) for p in all_files_0],
    key=lambda slice: slice.SliceLocation,
)

fig = plt.figure(
    figsize=(figure_size[0], figure_size[1]), constrained_layout=True)

gs = gridspec.GridSpec(
    nrows=2, ncols=6, figure=fig, wspace=0.08, hspace=0.04, height_ratios=[1, 1]
)

ax_1 = plt.subplot(gs[0, :2])
ax_2 = plt.subplot(gs[0, 2:4])
ax_3 = plt.subplot(gs[0, 4:])
ax_4 = plt.subplot(gs[1, 1:3])
ax_5 = plt.subplot(gs[1, 3:5])
axs = [ax_1, ax_2, ax_3, ax_4, ax_5]

for slice_counter, ax in enumerate(axs):
    DICOM_2D_image(
        ax,
        "bone",
        mri_data_0[slice_counter],
        f"MRI {get_plane(mri_data_0[slice_counter])} plane image at "
        f"{mri_data_0[slice_counter].SliceLocation:^.1f} mm slice location",
        fontsize_adjustment=-1,
        patien_age="52",
    )

fig.suptitle(
    "Visualization Comparison of DICOM Images in Different Planes",
    fontsize="x-large",
    x=0.5,
    y=-0.02,
)

plt.show()

NIfTI in Python¶

Extraction of NIfTI File¶

In [26]:
path_to_dicom = (
    "../Datasets/An MRI Dicom Data Set of the Head of a Normal Male Human Aged 52/"
    "ST000000/SE000001/"
)
output_folder = (
    "../Datasets/An MRI Dicom Data Set of the Head of a Normal Male Human Aged 52/"
    "ST000000/"
)
# The `dicom2nifti.convert_dir` module provides the `convert_directory` function, which is
# automatically imported through the package initialization, to sort all DICOM files one-by-one
# by series, and then convert them to NIfTI format

# By default, the `convert_directory` function is set to True for the `compression` parameter,
# which is used to enable or disable Gzip compression, and for the `reorient` parameter,
# which is used to reorient the DICOM images according to the LAS direction
dicom2nifti.convert_directory(path_to_dicom, output_folder)

# NIfTI stands for Neuroimaging Informatics Technology Initiative and is co-sponsored by
# the National Institute of Mental Health and the National Institute of Neurological Disorders
# and Stroke

# NIfTI defines a neuroimaging data file format designed to meet the needs of the fMRI research
# community, and although some software continues to use other file formats, NIfTI has become
# the most widely used storage standard for fMRI and other MRI research data files

# NIfTI images can be compressed using a standard open-source algorithm called Gzip, which
# greatly reduces file size and therefore the amount of storage required for imaging data

# NIfTI defines a single image file ending in the `.nii` extension, while the compressed
# version appends the `.gz` extension, i.e. `.nii.gz`
nifti = nib.load(f"{output_folder}201_t2w_tse.nii.gz")

tabulation = Form_Generator()
tabulation.heading_printer("Getting and Extracting NIfTI Files by Conversion")

statements = [
    """
path_to_dicom = (
    "../Datasets/An MRI Dicom Data Set of the Head of a Normal Male Human Aged 52/"
    "ST000000/SE000001/"
)
output_folder = (
    "../Datasets/An MRI Dicom Data Set of the Head of a Normal Male Human Aged 52/"
    "ST000000/"
)
dicom2nifti.convert_directory(path_to_dicom, output_folder)

nifti = nib.load(f"{output_folder}201_t2w_tse.nii.gz")
"""
]
tabulation.statement_generator(statements)

variables = ["nifti"]
values = [str(nifti)]
tabulation.variable_generator(variables, values)

expressions = [
    'nifti.header["qoffset_x"]',
    "nifti.header.get_data_shape()",
    "nifti.shape",
]
results = [
    str(nifti.header["qoffset_x"]),
    str(nifti.header.get_data_shape()),
    str(nifti.shape),
]
tabulation.expression_generator(expressions, results)
Getting and Extracting NIfTI Files by Conversion

    +-------------------------------------------------------------+
    | Statement                                                   |
    +-------------------------------------------------------------+
    | path_to_dicom = (                                           |
    |     "../Datasets/An MRI Dicom Data Set of the Head of a     |
    |     Normal Male Human Aged 52/"                             |
    |     "ST000000/SE000001/"                                    |
    | )                                                           |
    | output_folder = (                                           |
    |     "../Datasets/An MRI Dicom Data Set of the Head of a     |
    |     Normal Male Human Aged 52/"                             |
    |     "ST000000/"                                             |
    | )                                                           |
    | dicom2nifti.convert_directory(path_to_dicom, output_folder) |
    |                                                             |
    | nifti = nib.load(f"{output_folder}201_t2w_tse.nii.gz")      |
    +-------------------------------------------------------------+
    +----------+--------------------------------------------------+
    | Variable | Value                                            |
    +----------+--------------------------------------------------+
    | nifti    | ⟨class 'nibabel.nifti1.Nifti1Image'⟩             |
    |          | data shape (256, 256, 27)                        |
    |          | affine:                                          |
    |          | [[-9.36898887e-01  3.20514254e-02                |
    |          |         6.37919828e-02  1.15272324e+02]          |
    |          |  [ 3.03588901e-02  9.27917123e-01                |
    |          |         -8.33337128e-01 -9.72956390e+01]         |
    |          |  [ 1.43172191e-02  1.29802674e-01                |
    |          |         5.94150448e+00 -8.23735046e+01]          |
    |          |  [ 0.00000000e+00  0.00000000e+00                |
    |          |         0.00000000e+00  1.00000000e+00]]         |
    |          | metadata:                                        |
    |          | ⟨class 'nibabel.nifti1.Nifti1Header'⟩ object,    |
    |          |         endian='<'                               |
    |          | sizeof_hdr      : 348                            |
    |          | data_type       : b''                            |
    |          | db_name         : b''                            |
    |          | extents         : 0                              |
    |          | session_error   : 0                              |
    |          | regular         : b''                            |
    |          | dim_info        : 0                              |
    |          | dim             : [  3 256 256  27   1   1   1   |
    |          |         1]                                       |
    |          | intent_p1       : 0.0                            |
    |          | intent_p2       : 0.0                            |
    |          | intent_p3       : 0.0                            |
    |          | intent_code     : none                           |
    |          | datatype        : uint16                         |
    |          | bitpix          : 16                             |
    |          | slice_start     : 0                              |
    |          | pixdim          : [-1.         0.9375     0.9375 |
    |          |             5.9999995  1.         1.             |
    |          |   1.         1.       ]                          |
    |          | vox_offset      : 0.0                            |
    |          | scl_slope       : nan                            |
    |          | scl_inter       : nan                            |
    |          | slice_end       : 0                              |
    |          | slice_code      : unknown                        |
    |          | xyzt_units      : 2                              |
    |          | cal_max         : 0.0                            |
    |          | cal_min         : 0.0                            |
    |          | slice_duration  : 0.0                            |
    |          | toffset         : 0.0                            |
    |          | glmax           : 0                              |
    |          | glmin           : 0                              |
    |          | descrip         : b''                            |
    |          | aux_file        : b''                            |
    |          | qform_code      : unknown                        |
    |          | sform_code      : aligned                        |
    |          | quatern_b       : -0.016685799                   |
    |          | quatern_c       : -0.9974202                     |
    |          | quatern_d       : -0.069515765                   |
    |          | qoffset_x       : 115.27232                      |
    |          | qoffset_y       : -97.29564                      |
    |          | qoffset_z       : -82.373505                     |
    |          | srow_x          : [-9.3689889e-01  3.2051425e-02 |
    |          |          6.3791983e-02  1.1527232e+02]           |
    |          | srow_y          : [ 3.035889e-02  9.279171e-01   |
    |          |         -8.333371e-01 -9.729564e+01]             |
    |          | srow_z          : [ 1.4317219e-02  1.2980267e-01 |
    |          |          5.9415045e+00 -8.2373505e+01]           |
    |          | intent_name     : b''                            |
    |          | magic           : b'n+1'                         |
    +----------+--------------------------------------------------+
    +-------------------------------+----------------+
    | Expression                    | Result         |
    +-------------------------------+----------------+
    | nifti.header["qoffset_x"]     | 115.27232      |
    | nifti.header.get_data_shape() | (256, 256, 27) |
    | nifti.shape                   | (256, 256, 27) |
    +-------------------------------+----------------+
In [27]:
# Similar to DICOM files, image pixel data from NIfTI files can be extracted using the
# `get_fdata` function
image_array = nifti.get_fdata()

tabulation = Form_Generator()
tabulation.heading_printer("Getting pixel data from NIfTI file")

statements = ["image_array = nifti.get_fdata()"]
tabulation.statement_generator(statements)

variables = ["image_array"]
values = [str(reprlib_rules.repr(image_array))]
tabulation.variable_generator(variables, values)

expressions = [
    "image_array.dtype",
    "image_array.shape",
]
results = [
    str(image_array.dtype),
    str(image_array.shape),
]
tabulation.expression_generator(expressions, results)
Getting pixel data from NIfTI file

    +---------------------------------+
    | Statement                       |
    +---------------------------------+
    | image_array = nifti.get_fdata() |
    +---------------------------------+
    +-------------+------------------------------------------+
    | Variable    | Value                                    |
    +-------------+------------------------------------------+
    | image_array | array([[[0., 0., 0., ..., 0., 0., 0.],   |
    |             |         [0., 0., 0., ..., 0., 0., 0.],   |
    |             |         [0., 0., 0., ..., 0., 0., 0.],   |
    |             |       ... ...,                           |
    |             |         [0., 0., 0., ..., 0., 0., 0.],   |
    |             |         [0., 0., 0., ..., 0., 0., 0.],   |
    |             |         [0., 0., 0., ..., 0., 0., 0.]]]) |
    +-------------+------------------------------------------+
    +-------------------+----------------+
    | Expression        | Result         |
    +-------------------+----------------+
    | image_array.dtype | float64        |
    | image_array.shape | (256, 256, 27) |
    +-------------------+----------------+
In [28]:
fig, axis = plt.subplots(3, 3, figsize=(10, 10))

slice_counter = 0
for i in range(3):
    for j in range(3):
        image_display(image_array[:, :, slice_counter], axis[i][j], " ",
                      cmap="bone")
        slice_counter += 1
In [29]:
from scipy import ndimage


def grid_NIfTI_2D_image(inputs, ax, n=None, row_size=9, pad_value=1, rotation=0):
    #
    inputs = inputs.transpose(2, 0, 1)
    if n is None:
        n = inputs.shape[0]
    nrows = math.ceil(n / row_size)
    cmap = "bone"

    rotated_inputs = inputs.copy()
    for i in range(len(rotated_inputs)):
        #
        rotated_inputs[i, :, :] = ndimage.rotate(
            rotated_inputs[i, :, :], rotation, reshape=True
        )
    inputs = torch.Tensor(rotated_inputs.astype(int)).unsqueeze(1)

    for i in range(nrows):
        if len(inputs) > row_size:
            row_images = inputs[:row_size]
            inputs = inputs[row_size:]
        else:
            row_images = inputs[:]
        grid_row = make_grid(
            row_images, nrow=row_size, normalize=True, pad_value=pad_value
        )
        scalar_row = grid_row.numpy().transpose(1, 2, 0)[:, :, 0]
        if i == 0:
            if nrows != 1:
                scalar_image = scalar_row
            else:
                scalar_image = np.full(
                    (
                        scalar_row.shape[0],
                        math.ceil(scalar_row.shape[1] / inputs.shape[0] * row_size),
                    ),
                    np.nan,
                )
                scalar_image[:, : scalar_row.shape[1]] = scalar_row
        elif scalar_image.shape[1] == scalar_row.shape[1]:
            scalar_image = np.concatenate((scalar_image, scalar_row), axis=0)
        else:
            image_to_concat = np.full(
                (scalar_image.shape[0] + scalar_row.shape[0], scalar_image.shape[1]),
                np.nan,
            )
            image_to_concat[: scalar_image.shape[0], :] = scalar_image
            image_to_concat[scalar_image.shape[0] :, : scalar_row.shape[1]] = scalar_row
            scalar_image = image_to_concat
    ax.imshow(scalar_image, cmap=cmap)
    ax.set(xticks=[], yticks=[], frame_on=False)
    ax.grid(False)
    return


fig, axs = plt.subplots(2, 1, figsize=(figure_size[0], figure_size[1] / 7 * 6))

grid_NIfTI_2D_image(image_array, axs[0], rotation=0)
grid_NIfTI_2D_image(image_array, axs[1], rotation=90)

fig.suptitle(
    "",
    fontsize="x-large",
    x=0.5,
    y=0,
)

plt.tight_layout()
plt.show()
- [Extraction of DICOM Files](#toc5_1_1_)
- [Visualization of 2D DICOM Files](#toc5_1_2_)
- [Visualization of 3D Data Stored as Multiple 2D DICOM Files](#toc5_1_3_)
- [Extraction of DICOM Files by Auto-Detection](#toc5_1_4_)